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Abstract 



We have investigated the structure of double quantum dots vertically cou- 
pled at zero magnetic field within local spin-density functional theory. The 
dots are identical and have a finite width, and the whole system is axially 
symmetric. We first discuss the effect of thickness on the addition spectrum 
of one single dot. Next, we describe the structure of coupled dots as a func- 
tion of the inter dot distance for different electron numbers. Addition spectra, 
Hund's rule and molecular-type configurations are discussed. It is shown that 
self-interaction corrections to the density functional results do not play a very 
important role in the calculated addition spectra. 
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I. INTRODUCTION 



The study of systems with a small number (N) of electrons confined to a quasi two- 
dimensional (2D) semiconductor quantum dot (QD) constitutes a subject of growing interest 
(see for instance Refs. [3]|| and references therein). One reason for this interest is that their 
electronic properties can be selected with some freedom by tailoring the shape of the lateral 
confining potential. In this sense, they are often referred to as artificial atoms. Recently, 
circular and elliptically disk-shaped QD's have been built in a very clean way, and their 
properties have been thoroughly studiedU function of N. 

Usually, QD's are described as true 2D systems. This seems justified, as it has been 
widely tested by comparing theory with experiment. However, this comparison may be 
obscured by the fact that the parameters defining the confining potential are adjusted and 
might mask 3D effects. However, complex 2D calculations have been carried out without 
further restrictions or imposed symmetriesJUil Full 3D calculations also exist in the Hartree 



or Hartree- Fock (HF) approximations for veryfew electrons (see for instance Refs. p~3| r p~T| ) - 
and also within the density functional theory.HMii2 

A systematic study of the role of dimensionality in QD structure has been presented 
by Rontani et al.SH^ In spite of the success of 2D models in describing the properties 
of QD's subjected to perpendicular magnetic fields, the width in the growth direction of 
most experimentally studied QD's is around one order of magnitude smaller than their 
typical radius. The effect of a small but finite vertical extension on the QD structure is 
worth studying considering that a recent extension of single QD studies to vertically coupled 
quantum dots" may render including the z extension of the constituent dots unavoidable 
to describe the experimental data.Bi 



or quantum dot molecules, have been the- 
Only in Ref. |2^ has the z extension of the 



Vertically coupled dots, also called artificii 
oretically addressed in a number of works.c2rE 
constituent QD's been taken into account; in the other references, it was neglected and con- 
sequently, their results cannot be reliable when the interdot distance is comparable to their 
z extension, which is an interesting physical situation.1'1 

In this work we present a study of vertically coupled, cylindrical QD's in the local 
spin-density functional theory (LSDFT). We restrict our description to axially symmetric 
configurations and identical QD's, and limit ourselves to zero magnetic field B. While these 
two latter conditions are easy to relax without increasing the amount of numerical work too 
much, axial symmetry breaking would require a more demanding 3D calculation. 

Our scheme is based on the application of the so-calledH imaginary-time method (ITM). 
It requires a discretization of the Kohn-Sham (KS) equations on a spatial mesh that can be 
straightforwardly implemented on a personal computer, and avoids expansion of the single- 
particle (sp) wave functions in a basis, large matrix diagonalizations and tests of the stability 
of the results against changes in the size of the basis. 

In contrast to standard methods for computing molecular structured, we do not postulate 
from the start that the artificial molecular orbitals can be expressed as linear combinations of 
artificial atomic orbitals (LCAAO) and consequently, we avoid making approximations such 
as complete neglect of differential overlap (CNDO) where the individual wave functions 
of electrons associated with different artificial atoms are taken as being orthonormal, as 
happens when the QD's are either two dimensional or far apart. 
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The present approach to the description of vertically coupled QD's constitutes a clear 
improvement on previous LSDFT calculations^ in which the electrons are located on one 
of the dots, and as a consequence, can only be electrostatically coupled, even when they lie 
at short distances. With respect to the generalized Hubbard model plus the diagonalization 
scheme of Ref. |22], the present improvement is the possible application of LSDFT to systems 
with large numbers of electrons. Moreover, one could think of a trivial generalization to 
non-identical dots, paying the well known token that LSDFT treats the exchange Coulomb 
term locally, and that LSDFT configurations are usually mixtures of many-electron states 
with the same value of the total spin projection S z , however, with different total spin S. For 
small systems, this drawback (which is also inherent to the generalized Hubbard anproac hi) 
can be removed with a shell-model calculation in a restricted sp valence space. 0^ Another 
difference between our approach and the Hubbard model is that even if the exchange energy 
is treated locally, it is not restricted to involving only electrons located on the same dot. 
This might be advantageous in the strong coupling case when the dots are close together 
and the quantum mechanical coupling due to electron exchange is important. 

This paper is organized as follows. In Section II we present the formalism and the 
essentials of the ITM. Results for one single thick QD and for two coupled dots are presented 
in Section III, and a summary is presented in Section IV. In an Appendix we present the self- 
interaction correction!! to the density functional results obtained for two extreme interdot 
distances and several electron numbers. 



II. METHODOLOGICAL APPROACH 

Within LSDFT, the ground-state of the system is obtained by solving the Kohn-Sham 
equations. The problem is simplified by the imposed axial symmetry around the z axis, 
which allows one to write the sp wave functions as (f) n i a (r, z,8,a) = u n i a {r, z)e~ lW Xa with 
/ = 0, ±1, ±2 . . ., where —I is the projection of the sp orbital angular momentum on the 
symmetry axis. 

We have used effective atomic units h = e 2 /e = m =1, where e is the dielectric constant, 
and m the electron effective mass. In units of the bare electron mass m e one has m = m*m e . 
In this system, the length unit is the effective Bohr radius = aoe/m*, and the energy unit 
is the effective Hartree H* = Hm*/e 2 . In the numerical applications we have considered 
GaAs, for which we have taken e = 12.4, and m* = 0.067. This yields ~ 97.94 A and 
H* ~ 11.86 meV. 

In cylindrical coordinates the KS equations read 




2 \dr 2 r dr r 2 dz 2 / 

(1) 

+ V H + V xc + W xc r] a ] u nla (r, z) = e nla u nla (r, z) , 

where 77^ = H- 1 ( — 1 ) for cr=|(|), V c f(r,z) is the confining potential, V H (r,z) is the direct 
Coulomb potential, and V xc = d£ xc (n,m)/dn\ gs and W xc = d£ xc (n,m)/dm\ gs are the vari- 
ations of the exchange-correlation energy density S xc (n,m) in terms of the electron density 
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n(r, z) and of the local spin magnetization m(r, z) = n^(r, z) — n^-(r, z) taken at the ground 
state (gs). 

As usual, £ xc (n, m) = S x (n, m) + S c (n, m) has been built from 3D homogeneous electron 
gas calculations. This yields a well-knownc3, simple analytical expression for the exchange 
contribution £ x (n,m). For the correlation contribution S c (n,m) we have used two different 
parameterizations, both based on the results of Ceperley and Alder.il The first was proposed 
by Vosko, Wilk and Nusairll, and the second by Perdew and Zunger.S We have checked 
that they yield the same results, and all results presented in this work have been obtained 
with the exchange-correlation energy density proposed by Perdew and Zunger. 

For a double QD the confining potential V c f(r,z) has been taken to be parabolic with 
frequency uq in the xy plane, plus a symmetric double quantum well in the z direction. For a 
single QD we have also used a parabolic confining potential in the xy plane, together with a 
quantum well in the z direction. Any other axially symmetric V c f(r, z) can be implemented 
as well. 

V H (r\z) was obtained solving the Poisson equation using the conjugate gradient 
methodcaliU (CGM). This requires a knowledge of V H (r,z) at the mesh boundary which 
can be obtained by direct integration. Due to axial symmetry 

fco z'+oo 

V H (r,z)=2 r'dr' dz' An(f ') [(r + r') 2 + (z - z') 2 ] 1/2 E{a 2 ) , (2) 

JO J-oo 

where E is the complete elliptic integral of the second kind0 and a 2 = 4rr'/[(r + r') 2 + (z — 

m 

We have discretized the KS and Poisson equations using fc-point Lagrange formulae for 
the r- and z-derivatives, and k + 1-point Lagrange formulae for the integrations.^ The mesh 
size has to be such that the discretized wave functions u(ri, Zj) at the mesh boundary can be 
safely taken as zero. This is one boundary condition for physically acceptable solutions to 
Eq. (D, the other one is the regularity of the sp wave functions at r = 0. In our scheme the 
Ar and Az steps may have different values. The high precision demanded by the calculation 
imposes restrictions on the possible values of k, Ar, and Az and will be discussed below. 
This space discretization scheme offers an efficient calculation of sp wave functions, thus of 
the electron densities and direct Coulomb and exchange-correlation potentials. 

The imaginary time method is described in detail in Ref. |28|. It is based on the observation 
that the discretized time-evolution operator 7i(t) for the time-dependent KS equations 



ut, — - 



which formally yields 



^ = «M (3) 



^) =exp(~Am (n) )4 n) , (4) 



<f +1) = exp {- % -AtH^\ 4 n) 

where At is the time step and n indicates the step iteration, for imaginary At = — iAt 
(Ar > 0) produces a decrease of the KS energy. In imaginary time, the wave functions 

are no longer orthonormal, and a Gram-Schmidt orthonormalization has to be 

(,(*»+!) l fL. nm r lT >+ 1 )- 



carried out after each iteration to obtain {m } from {\EQ }. To first order in Ar, Eq. 
(ft) becomes 
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^(x.zq^y, (5) 

which shows the simplicity of the method and its practical implementation: after discretizing 
the KS Hamiltonian and wave functions, it essentially reduces to repeated application of the 
KS Hamiltonian to the previous wave functions. Moreover, it is easy to check from Eq. @ 
that upon convergence 



e 



(n+l) = h /, C'i hT( i,>+li\ 

j At 



*f +1) ) (6) 



coincide with the KS sp energies. Equation @ shows that the ITM belongs to the general 
class of relaxation methods employed to solve partial differential equations. This provides 
a simple criterion for fixing At in such a way that the imaginary time evolution is stable^, 
namely h 2 AT max /(2mA 2 ) < 1/4, with A the smallest between the r and z steps. In actual 
calculations we have taken At = 0.1 AT max . 

To start the ITM iteration one needs a set of sp wave functions and energies to build the 
initial sp level scheme. We have used two such sets. The first consists of the wave functions 
of an axially symmetric, 3D harmonic oscillator potential of frequencies Uq and u z in the 
radial and z directions, respectively. This potential gives rise to analytical solutions even 
in the presence of a constant magnetic field in the z direction, and is a suitable confining 
potential for a single QD. In this case, the sp Hamiltonian is separable and the wave function 
can be written as u(r, z) = TZ(r) Z(z)/ 'y2~7r with 



Zlz) 




2 



2a 



(7) 



where a = W^/2mw , C = \j mijJ zl^i and and H nz are generalized Laguerre and Hermite 
polynomials^, respectively. The sp energies are £ nr ,n z ,i = huJ [2n r + |/| + 1] + huj z [n z + 1/2] 
with n r and n z equal to 0, 1, . . . 

For a double QD, we have found it convenient to choose Z(z) as the lowest energy 
eigenfunctions of a ID double quantum well. The QD thickness is such that for not too 
many electrons, only the two lowest states are needed. If the double dot is symmetric, these 
solutions are either even or odd under reflection z — > —z. Actually, for one QD, Z(z) can 
also be the wave function in the ID quantum well. 

Both for a single QD and for a symmetric double dot, the single-electron wave functions 
are characterized by the values of l z , s z and parity, i.e., they must be either symmetric or 
antisymmetric under inversion f — > — r, and be either even or odd when reflecting z — > —z. 
All these symmetries are included in the starting separable wave functions. Indeed, the 
Hermite polynomials are even or odd depending on n z , and the parity of a sp level is simply 
(_)n z +|*]_ Ag the K g Hamiltonian and the ITM preserve these symmetries, in the course 
of the iteration procedure the sp wave functions, which are no longer separable, keep their 
initial quantum numbers (orbital and spin projection on the z axis, parity and reflection 
symmetry with respect to the z = plane) which are conserved quantities. 
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To check the numerical scheme we have carried out extensive and systematic tests. A 
first test on the discretization and iteration procedure consisted in numerically solving an 
axially deformed harmonic oscillator potential. We have exactly reproduced the spectrum 
£n r ,n z ,i given after Eqs. (0). The implementation of the CGM has been successfully tested 
comparing results computed for a spherical Gaussian charge distribution with the analytical 
distributions. 

As another test of the numerical code, we have compared the total energy calculated 
from a straightforward integration of the energy density with the expression in terms of the 
sp energies derived from the KS equations. Writing the correlation energy densityil'i! as 
S c [n,m] = nScorln, £], where S cor [n,^] is the correlation energy per electron and £ = m/n is 
the local spin polarization, one obtains 

E = E 3 - / df | \ V H (f) n(r) + ± S x [n, m) + d£ ^® ^ } | (g) 

We have checked that in the worst case, the energies calculated with either method agree 
to within one part in 10 4 for our combination of k = 7 formulae and the values Ar = Az = 
0.12 <2q used in this work.S Simpler three point formulae (k = 3) turned out to be inaccurate 
for reasonable spatial steps, and k = 9 and 11 formulae, which allow larger steps, did not 
appreciably reduce the computing time. It might be interesting to remark that most of this 
disagreement arises from the integral of the external confining potential for a quantum well 
in the z direction, due to the sharp discontinuity in the given potential. We have checked 
that for a parabolic confining potential in the z direction, both methods of evaluating the 
total energy agree to within one part in 10 7 . 

In the case of double QD's, we have also checked that the results coincide irrespective 
of whether we start the iteration from pure harmonic oscillator wave functions, or from 
the better choice of double quantum well eigenfunctions in the z direction. For even- A 
systems we have always started the iteration with non-identical sp potentials for spin-up 
and spin-down electrons, to avoid artificial configurations with S z = 0. 

III. RESULTS 

A. Single quantum dot 

We have first addressed the addition energies of one quantum dot hosting up to N = 21 
electrons. The confining potential in the z direction is a quantum well W = 12 nm wide, 
which corresponds to the experimental well@, and Vq = 200 meV deep.0 Electrons are 
laterally confined by the parabolic potential rraJ^r 2 /2 for which we have tried different u 
values. 

Figure [1| shows the addition spectrum AA(N) = E(N + 1) - 2E(N) + E(N - 1) for u 
= 3, 5, and 10 meV. In addition to local maxima at shell filling values N = 2, 6, 12, and 20, 
other peaks appear at half-filling values N = 4 and 9. This is a consequence of Hund's rule, 
which establishes that degenerate electronic states in a shell are filled with parallel spins up 
to half-shell, in order to maximize the exchange interaction. 

Up to N = 12 the results look very similar to the experimental values!, especially for 
medium and weak r-confinement. For larger N values we have found a conspicuous even-odd 



6 



effect, instead of the weak maxima experimentally found at iV = 20, and especially at iV 
= 16. This result does not mean that Hund's rule is violated within this shell. Indeed, we 
have found spin alignment up to iV = 16, but the associated energy gain is not enough to 
produce the local maximum at N = 16. It is worth pointing out that the energies involved 
in the definition of AA(N) are very large, compared with the second energy difference. For 
example, for uq = 5 meV, E(16) = 1.048 eV, whereas AA(16) is ~ 3 meV. 

For strong r-confinement, our results are quite similar to those in Ref. [19], except for 



medium and small values. This is not surprising, since in this Ref. the 3D electron-electron 
energy is treated in Brst-crder perturbation theorjH, thus one should not expect it to hold 
if the confinement is not strong enough. The two-dimensional E(N + 1) — E(N) results at 
B — coincide with those in Ref. ^ when uj® = 5 meV. 

Our 3D dot is rather strongly confined in the z direction, and the results for AA(N) are 
similar to those obtained for pure 2D dots. We have computed the addition energies using 
the 2D dot model described in detail in Ref. ^D|; as shown in Fig. 2D and 3D results 
are very similar and qualitatively better than those from other different 2D models. As 
one increases the confinement in the z direction up to Vq = 300 meV, the addition energies 
become indistinguishable in the scale of Fig. [I]. We shall see that this is not the case for 
double dots, and that even qualitative differences appear if the dots are strongly coupled. 

In spite of this quasi two-dimensional behaviour, the electron density spreads beyond the 
well up to distances that are relevant for vertically coupled dots. This is illustrated in Fig. 
[3], where we have plotted the n(z) electron density defined as 



CO 



n(z) = / drrn(r,z) (9) 



o 



corresponding to N= 12, W— 12 nm, and Vq = 200 and 300 meV [notice that iV = 
27r / dzn(z)}. The electron densities spill out of the quantum well, the smaller the con- 
finement, the larger the effect. 

Finally, in Fig. ^ we represent the density of a 2D dot hosting 12 electrons, and the 
n(r) density of the corresponding 3D dot defined as n(r) = J dzn(r,z) [notice that N = 
2-7T / drr n(r)], both for ujq = 5 meV. These radial densities are very similar. The 3D densities 
n(r) for Vq = 300 meV, instead of 200 meV, are indistinguishable within the scale of the 
figure. 

These comparisons allow one to infer that the experimental dots are quasi two- 
dimensional systems to a large extent, with moderate lateral confinement, ujq < 5 meV. 



B. Double quantum dots 

We have modelled a symmetric double dot by a parabolic confining potential with fre- 
quency ujq = 5 meV in the r direction, and a symmetric double quantum well in the z 
direction. Each quantum well has a width equal to 12 nm, and is separated from the other 
by a barrier of thickness d that varies from 1 to 9 nm. Some experimental results are 
available! for this system at d > 2.5 nm. Results for well depths Vq = 200 and 300 meV will 
be discussed. 

Figure [| shows the n(z) density profiles corresponding to a quantum- mechanically cou- 
pled configuration (d = 2.5 nm), and to an electrostatically coupled configuration (d = 7.5 



7 



nm). Only in the former case, are the electrons delocalized. A quantitative measure of 
this localization is provided by the energy splitting between symmetric and antisymmetric 
sp states (more precisely between even and odd states with respect to specular reflection 
z — > —z), As as- Indeed, for N = 1 and large distances, the symmetric and antisymmetric 
states are degenerate and Asas approaches zero. As d decreases, the coupling increases and 
so does Asas- This effect depends weakly on N. We have plotted Asas as a function of d 
in Fig. for the lowest I = sp levels of systems with N = 1 and 20. 

Due to the characteristics of our double quantum well, Asas depends exponentially on d, 
AsAs(d) = A exp(-d/d ). For V = 200 meV we have (d , A„) = (1.79 nm, 19.2 meV) for 
N — 1 and (1.76 nm, 17.6 meV) for N = 20, whereas for V = 300 meV, (d , A ) = (1.44 nm, 
17.9 meV) for N = 1 and (1.42 nm, 16.5 meV) for N = 20. These values are similar to those 
in Refs. ^ and ^71 It can be seen from Fig. |5] that for a given interdot distance, enlarging Vq 
decreases the coupling between the dots, as the density overlap diminishes. Figure [7| shows 
the densities n(r, z) for N = 12 and for two d values corresponding to quantal (d = 2.5 nm) 
and to electrostatic (d = 7.5 nm) coupling. 

Figure H shows the addition spectrum for quantum mechanically coupled (d = 2.5 nm), 
and electrostatically coupled (d = 7.5 nm) dots, as well as the results corresponding to 
the single 3D dot. They have been obtained for Vq = 200 meV. Figure |] shows the same 
spectrum for V = 300 meV. As expected, changes mostly appear in the strong coupling 
case (d = 2.5 nm). Unexpectedly enough, these changes are qualitative, with, for example, 
maxima in AA(N) changed into minima. 

Fig. H indicates that at d = 7.5 nm the dots are well apart, only influenced by the 
electrostatic coupling, and the addition energies are insensitive to the value of Vq. There 
are no experimental results in the literature for AA(N) at this interdot distance. However, 
the experimental analysis^ of the derivative of the drain intensity with respect to the drain 
voltage versus drain voltage indicates that the electrons in the dots are indeed delocalized 
for d = 2.5 nm, and rather localized for d = 7.5 nm. 

The weak coupling at d — 7.5 nm allows us to interpret the appearance of several peaks 
in AA(N), such as those at N = 4, 8 and 12, as due to the fact that 2, 4, and 6 electrons 
on each dot already yield maxima in the addition spectrum of a single QD. The remaining 
peak at N = 2 corresponds now to half-filling the first shell in each dot, which would close 
at iV = 2. This is caused by the localization of one electron on each constituent dot.0 

Experimental results have been published!'! for d = 2.5 nm. We have not attempted to 
use Vq as a fitting parameter, so when comparing with experiment one should bear in mind 
both the sensitivity of AA(N) on the value of Vq in the strong coupling limit and the results 
displayed in Figs. R and ^. Published calculations correspond to depths lying in between 
these two values. IHHEZI 

One can see that maxima of AA(N) decrease on the average when compared to the 
isolated QD. This is in agreement with experiment. Notice however that the experimental 
AA(N) for the strongly coupled double dotu is approximately half that corresponding to the 
single dot, especially when N > 10. This feature is not reproduced by the calculations. 

Globally, the results for d = 2.5 nm are better reproduced with Vq = 200 meV up to 
N = 12, and with Vq = 300 meV for larger N values. Probably, a value of Vq in between 
might improve the agreement with experiment. Other possibilities such as an iV-dependent 
uoq might also be consideredi@, or even some asymmetry in the double well.B. We have not 
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tried these possibilities^, but we have checked that self-interaction corrections (SIC), which 
are usually not included in these kind of calculations, do not change the addition spectrum. 
The results are presented in the Appendix. 

For a given electron number, the gs configuration may change as a function of the barrier 
thickness. The new 'phases\ i.e. gs configurations which appear as a function of d have 
been thoroughly discussedMS To label them, we have adopted the standard convention 
of molecular physics for sp electronic orbitals as a, 7r, 6, . . . , if I = 0, 1, 2, . . ., and upper 
case Greek letters are used for the total orbital angular momentum. We have also used an 
adapted version!! of ordinary spectroscopic notation 2S+1 L^ U where S is the total \S Z \, and 
L is the total \L Z \. The superscript +(— ) refers to even (odd) states under reflection with 
respect to the z = plane, and the subscript g{u) refers to positive (negative) parity states. 

We show the evolution with barrier thickness of the energy and gs molecular configuration 
for several N values in Fig. [Hj The vertical lines have been drawn to guide the eye, and 
different symbols have been used to identify different phases. All panels in the figure display 
some common trends. Initially, E(d) increases with d. The reason is twofold. On the one 
hand, at small d all occupied sp levels are specularly symmetric about the z = plane ('+' 
states), the specularly antisymmetric levels ('— ' states) lie at much higher energies (see Fig. 

On the other hand, the energies of the symmetric and antisymmetric state respectively 
increase and decrease with d, and eventually both states become degenerate at large interdot 
distances. This is a well-known feature of the one particle, one-dimensional double quantum 
well problem (see insert in Fig. §), which remains valid in the interacting many-electron 
calculation. The first phase transition takes place when the first antisymmetric sp state 
becomes occupied. At large distances, E(d) slowly decreases with d due to the decrease of 
the interdot Coulomb energy. These trends are also present in the Hubbard-like calculations 
of Ref. but only the lowering of E(d) due to the interdot Coulomb energy is qualitatively 
reproduced by the calculations of Ref. |27|, which fail to yield the energy growth at short 
distances. 

In spite of the difference between the values of the gs energies reported in Refs. |22|j27| , 
and also with respect to the present work, which has to be mostly attributed to the different 
confining potentials in each calculation, up to N = 6 the phases are the same but appear 
at different d values. In this respect, our results are in closer agreement with those of Ref. 
p7| , possibly because the radial frequencies u are similar in both calculations. We wish to 
point out that the phase diagram obtained for Vq = 300 meV is qualitatively similar to the 
one shown in Fig. |TD|, but the phase transitions are shifted ~ 0.5 nm to the left. 



IV. SUMMARY 

We have used local spin-density functional theory to investigate the zero magnetic field 
structure of one single and two identical, vertically coupled QD's of finite thickness. While 
for one single dot, whose thickness corresponds to that of actual experimental devices, the 
addition spectrum is quite similar to that predicted by purely two-dimensional models, in 
the case of double dots, their vertical extension, is essential for a quantitative description of 
their quantum coupling, which influences the addition spectrum at short distances. For one 
single dot the calculated addition spectrum compares well with experiment, whereas for two 
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coupled dots the agreement is qualitative. This possibly reflects the fact that in the latter 
case the spectrum is more sensitive to the actual form of the bare confining potential. 

The phase sequence of ground state configurations which appear as a function of interdot 
distance is quite similar to that found in previous works§!izl, evolving from the 'atomic phase' 
of two strongly coupled dots to the 'atomic phase' of two weakly coupled dots through a 
series of 'molecular-type phases' at intermediate distances. This is a rather robust picture, 
as it arises from the underlying single-electron structure of the bare confining potential. 
Indeed, the vertical confinement is so strong that at short distances only symmetric sp 
states are occupied, the antisymmetric ones lying at quite high energies. This originates 
the atomic phase of two strongly coupled dots. As the interdot distance increases the 
Asas g a P decreases, and the symmetric and antisymmetric sp states eventually become 
degenerate, originating the atomic phase of two weakly coupled dots. The molecular-type 
configurations appear at intermediate distances where A$as is similar to the other energy 
scale of the system, namely hu , and when the number of electrons is large enough so 
that the system can minimize its energy populating antisymmetric states. The larger the 
number of electrons, the larger the number of populated antisymmetric states. This causes 
the number of intermediate molecular phases to increase with N. 

In spite of the mentioned qualitative agreement with previous results, the calculated 
phase diagrams are as sensitive to the shape of the bare confining potential as the addition 
spectra. It would be desirable that any prediction of the actual appearance of the phase 
diagrams should be based on a model that describes, at least qualitatively, the corresponding 
experimental addition spectrum. 
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APPENDIX A: 

It is well known that the 'exact' density functional for the gs energy is self-interaction-free 
(see for instance Ref. HBD, but it is not the case of its current approximations, such as LSDFT. 



One possible way of removing this drawback is to use the self-interaction correction (SIC) 
proposed by Perdew and Zungeiil, which introduces an orbital-dependent single-particle 
potential, that improves the total energy of the electronic system and yields sp eigenvalues 
which approximate the physical removal energies more closely, at the price of rendering the 
KS minimization far more cumbersome. Since SIC's are relatively more important for few 
electron systems, which is the present case, we have tested their effect on the results shown 



in the body of this paper, in two extreme configurations. We refer the reader to Refs. B2 



and 46 for a thorough description of SIC. Here we only give the essential details of the 



application of the method to our physical problem. 
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Within the method of Perdew and Zunger, the sp potential in Eq. ([]]) becomes orbital- 
dependent with the change 



V H + V xc + W xc r] a — > 

V eff = V H + V xc + W xc r] a 

- V H [n nla \ - V xc [n nla , n nla \ - W xc [n nla , n nla ] r\ a , (Al) 

where n nla = \u n i a (r, z)\ 2 is the 'orbital density' and V H [n n i a ] is obtained solving the Poisson 
equation &.V H [n n i a ] = —4tt rin^^r, z) as indicated in Section II. After selfconsistency is 
achieved, the total energy can be obtained from Eq. (|8|), which has been written as 

E — > E - ^ $Z J drV H [n n i a ] n n i a {r ) - ^ J dr £ xc [n n i a , n n i a ] . (A2) 

nlcr nla 

We have drawn in Fig. |ll] the difference between SIC and LSDFT addition energies 
defined as n(N) = E(N) — E(N — 1) corresponding to the Vq = 200 meV double quantum 
dot. The average difference is small, of the order of 0.4-0.5 meV. This difference almost 



cancels out in the addition spectrum AA(N), as can be seen in Fig. |T|. It constitutes an 
interesting result in itself, indicating that if one only wishes to obtain the addition spectrum, 
LSDFT does not need to be corrected for self-interaction effects. 
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FIGURES 



FIG. 1. Addition energies A A as functions of A" for a 3D single dot with Vq = 200 meV, W 
= 12 nm and different cuq values. 

FIG. 2. Addition spectrum as a function of A" for a 2D single dot and different ojq values. 

FIG. 3. n(z) densities [K)" 1 ] for a 3D single dot with N = 12, width W = 12 nm and well 
depths Vq = 200 and 300 meV. The vertical lines indicate the limits of the quantum well. 

FIG. 4. n(r) densities [(oq)" 2 ] for 2D and 3D single dots with N = 12 and ujq = 5 meV. The 
3D well width and depth are W = 12 nm and Vo = 200 meV respectively. 

FIG. 5. n{z) densities [(a*^ 1 } for an JV = 12 double dot of width W = 12 nm, well depth 
Vq, and barrier thickness d. Top panel, d = 7.5 nm. Bottom panel, d = 2.5 nm. The vertical lines 
indicate the limits of the double quantum well. 

FIG. 6. AsAs(d) for the lowest I = sp levels of a double dot with N = 1 and 20. The insert 
shows the energies E$ and Eas defining As as, As as = Eas — Es, for one of the cases presented 
here. 

FIG. 7. n(r, z) densities [(ao) -3 ] for a double dot with N = 12, width W = 12 nm, well depth 
Vq = 200 meV, and barrier thickness d. Top panel, d = 2.5 nm. Bottom panel, d = 7.5 nm. 

FIG. 8. Addition spectrum as a function of for a 3D single dot (circles), and for two coupled 
dots at d = 2.5 nm (squares) and d = 7.5 nm (diamonds). The depth of the double quantum well 
is V = 200 meV. 

FIG. 9. Same as Fig. | for Vq = 300 meV. 

FIG. 10. Energy and ground state molecular configurations of the double dot as a function of 
the barrier thickness for A^ = 3 to 7. The depth of the double quantum well is Vo = 200 meV. 

FIG. 11. Difference between SIC and LSDFT addition energies n{N) = E(N) - E(N - 1) for 
a double quantum dot of Vq =200 meV depth and barrier thickness d for N = 1 to 13. 

FIG. 12. Addition spectrum as a function of A^ for two coupled dots at d = 2.5 nm (squares) 
and d = 7.5 nm (diamonds). The depth of the double quantum well is Vq = 200 meV. Solid symbols 
represent LSDFT results, and empty symbols the results obtained after SIC's have been included. 
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